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EXECUTIVE  SUMMARY 


This  report  describes  the  formulation  of  NRL’s  hydrodynamic  (isopycnal)  nonlinear,  primitive 
equation,  ocean  circulation  model  in  a  spherical  layer  suitable  for  a  large  basin  or  the  global  domain. 
This  document  should  be  read  in  conjunction  with  NOARL  Report  35,  December  1991,  The  Navy 
Layered  Ocean  Model  User's  Guide  (hereafter  NOARLS5)  [25]. 

Spherical  geometry  requires  careful  treatment  of  the  equations  of  fluid  dynamics  to  ensure 
spurious  terms  are  not  introduced  during  the  layer  averaging  process,  particularly  in  the  momentum 
diffusion  expressions. 

Prom  these  layer  average  equations,  flnite  difference  approximations  are  derived  on  a  grid  of 
uniform  intervals  in  longitude  and  latitude  (not  necessarily  the  same).  The  resulting  flnite  differ¬ 
ence  equations  may  be  solved  eflBciently  on  modern  axchitectmre  scientific  computers  to  describe 
acciuately  ocean  cmrents,  temperature  and  salinity  distributions  over  large  time  and  space  scales. 
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Formulation  of  the  NRL  Layered  Ocean  Model 
in  Spherical  Coordinates 


1.  PRELIMINARIES 

1.1  Motivation 

Ocean  Modeling  has  evolved  in  20  years  from  simple  limited  resolution,  limited  area,  ideal 
domain,  Cartesian  geometry  computations  [7],  [21],  [8]  to  realistic  basin  scale  eddy  resolving  ocean 
models  [6],  [10],  [9].  This  progress  has  reflected  both  the  growth  in  computing  capacity  and  in 
the  sophistication  and  experience  of  numerical  oceanographers.  However,  the  actual  models  cannot 
form  a  continuous  range  as  the  transition  from  Cartesian  to  spherical  geometry  results  in  substantial 
changes  in  the  appearance  of  the  vector  calculus  operators  and  of  the  equations  of  fluid  flow,  in 
spite  of  the  coordinate  invariance  of  the  underlying  physics  and  mathematics.  This  change  in 
appearance  reflects  the  difficulties  of  expressing  the  laws  of  physics  in  a  curvilinear  coordinate 
system,  particularly  the  conservation  of  vector  quantities  such  as  momentum. 

Another  difficulty  is  that  while  the  choice  and  orientation  of  {x,y,z}  to  label  the  orthogonal 
Cartesian  coordinates  appears  to  be  common  across  most  physical  disciplines,  no  such  consensus 
exists  for  the  choice  and  orientation  of  spherical  or  spheroidal  coordinates.  While  most  mathemat¬ 
ical  and  theoretical  studies  favour  co-latitude  as  the  declination  coordinate,  many  oceanographers 
and  meteorologist  prefer  latitude.  However,  introductory  mathematics  and  physics  text  books 
tend  to  opt  for  the  former  choice.  This  can  lead  to  confusion  and  error.  The  labels  assigned  to  the 
three  spherical  surface  coordinates  also  vary  widely,  even  within  the  meteorology  and  oceanography 
literature. 

For  simplicity,  the  criteria  of  minimizing  diflFerences  between  the  Cartesian  and  spherical  models 
will  guide  the  choices  of  notation  and  orientation  used  in  this  presentation. 

A  note  of  caution.  Care  will  need  to  be  taken  to  distinguish  two-dimensional  mathematical 
operations  defined  on  the  surface  of  a  sphere  and  three-dimensional  physical  and  mathematical 
processes  defined  within  a  thin  spherical  layer.  The  purpose  of  the  Shallow  Water  set  of  approxi¬ 
mations  is  to  reduce  the  equations  that  apply  to  three-dimensional  physical  flow  throughout  a  thin 
fluid  layer  to  equations  describing  a  two-dimensional  model  flow  on  a  spherical  surface  in  a  way 
that  preserves  as  much  of  the  physics  of  the  3-D  flow  as  possible. 

1.2  Coordinates  and  unit  vectors. 

To  conform  to  the  conventional  right-handed  (rr,  y,  z)  Cartesian  coordinate  system  locally,  we 
chose  longitude,  latitude  and  radius  as  our  surface  spherical  coordinate  directions  in  that  order. 
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We  label  the  longitude  and  latitude  angles  and  6  respectively.  We  assume  that  0  increases  in 
am  eastward  direction,  and  that  0  increases  in  a  northward  direction  from  a  value  of  zero  at  the 
equator,  taking  negative  values  in  the  Southern  hemisphere.  The  two  spherical  surface  coordinate 
singularities  (the  North  and  South  Poles)  are  assumed  to  coincide  with  the  rotation  axis  of  the  earth, 
unless  stated  otherwise.  The  origin  of  the  longitude  [4>)  coordinate  can  be  chosen  for  convenience, 
with  the  Greenwich  meridian  assumed  to  be  the  default.  The  surface  of  the  earth  is  chosen  to  he 
on  a  sphere  at  a  constant  radius  a.  The  radieil  coordinate  is  chosen  to  increase  outwards.  We  will 
label  this  coordinate  r.  As  the  angles  (f>  and  0  are  dimensionless,  physical  distances  will  be  given 
implicitly  in  the  form  a<f>  or  a0  ,  where  o  is  a  suitable  mean  reference  radius  for  the  earth. 

We  label  the  unit  vectors  that  are  locally  tangent  to  the  surface  and  parallel  to  lines  of  constant 
latitude  and  constant  longitude  as  and  Sg  respectively.  We  assume  that  points  eastward 
and  €$  points  northward.  There  is  a  local  unit  normal  vector  Cr  pointing  out  of  the  sphere  at  each 
point  {(f>,  0,  a)  on  the  spherical  surface. 

The  distinction  between  a  sphere  and  true  geopotential  surfaces  will  be  ignored.  Gill  [4]  pp  91  - 
92  contains  a  full  discussion  of  the  error  incurred  in  approximating  the  shape  of  the  earth  by  a 
sphere.  There  he  estimates  that  a  maximum  error  of  0.17%  occiurs  in  the  geometrical  factors  in  the 
various  equations  if  the  earth’s  figure  is  approximated  as  a  sphere  instead  of  an  oblate  spheroid. 
Thus,  although  in  reality  horizontal  fluid  flows  are  locally  parallel  to  geopotential  surfaces,  very 
little  error  is  made  approximating  these  complicated  surfaces  by  concentric  spheres. 

1.3  History 

The  derivation  of  the  sphericad  Shallow  Water  equations  from  the  full  3-D  Navier  Stokes  equa¬ 
tions  in  a  spherical  coordinate  system  has  been  carried  out  by  several  authors  in  the  past  [11], 
[17].  While  early  agreement  was  reached  as  to  the  proper  form  for  the  momentum  advection  and 
conservation  of  mass  expressions,  opinions  continue  to  difier  as  to  the  correct  form  for  the  layer 
averaged  viscous  dissipation  in  the  presence  of  layer  thickness  and  horizontal  density  variations. 
We  followed  Wajsowicz  [24]  and  Shchepetkin  &  O’Brien  [18]  in  the  derivation  of  these  terms  and 
can  demonstrate  that  our  forms  generate  no  spurious  forces  for  analytic  test  cases  such  as  rigid 
rotation  of  a  varying  thickness  shell. 

Given  the  complexity  of  the  appearance  of  these  equations  in  spherical  coordinates,  many 
mathematically  equivalent  expressions  are  possible.  These  variations  in  appearance  can  lead  to 
variations  in  the  finite  difference  approximations  used  when  replacing  these  partial  differential 
equations  by  difference  equations.  These  can  give  rise  to  widely  differing  answers  to  the  same 
imderlying  mathematical  problem.  We  note  that  while  the  equations  given  in  Section  2  are  mathe¬ 
matically  equivalent  to  those  appearing  in  Shchepetkin  and  O’Brien  (1996)  [18],  the  most  recently 
published  contribution  to  this  discussion  that  we  are  aware  of;  they  are  not  identical.  They  differ 
in  the  form  taken  for  the  viscous  dissipation,  although  they  adhere  to  the  requirement  that  this 
term  be  the  divergence  of  a  symmetric  tensor.  The  form  we  use  is  taken  from  Wajsowicz  (1994) 
[24]. 
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2.  SPHERICAL  SHELL  FLUID  DYNAMICS 

The  reduction  of  the  three  dimensional  partial  differential  equations  describing  the  conservation 
of  mass,  momentum,  energy  and  matter  (tracers)  at  each  point  in  am  ocean  layer  to  a  set  of  two- 
dimensionaJ  partial  differential  equations  for  flow  on  a  spherical  layer  is  done  by  ntegrating  these 
equations  in  the  radial  direction  over  the  thickness  of  the  spherical  shell,  h  from  r  =  a  tor  =  a  +  h. 
This  averaging  process  requires  application  of  boundary  conditions  at  the  top  amd  bottom  of  eaich 
layer  and  the  conversion  of  surface  stresses  into  body  forces.  In  our  equations  these  are  subsumed 
into  the  “sources  —  sinks  ”  terms. 

In  Cartesian  coordinates,  the  integration  is  straightforward  and  yields  the  equations  set  out 
in  NOARLS5  [25]  or  Hurlbmt  amd  Thompson  [8].  On  a  spherical  shell  the  averaging  must  be 
done  with  care  to  avoid  generating  spurious  sources  or  sinks  of  energy,  mass  or  momentum  from 
curvature  terms.  After  a  lively  discussion  in  the  literatmre  (sf.  [1],  [2],  [5],  [11],  [12],  [13],  [15]  [17] 
[20]  [23]  &  [26])  a  consensus  has  emerged  as  to  the  correct  forms  for  these  equations  amd  we  have 
chosen  a  representation  that  closely  resembles  the  one  used  in  our  earlier  Cartesiam  work  [25]. 

2.1  2-D  Sphericad  Shell  Velocity  Field 


Classically,  the  derivation  of  the  Shallow  Water  equations  assumes  each  layer  is  homogeneous 
in  the  vertical  direction  amd  that  vertical  velocities  within  the  layer  cam  be  ignored  except  to 
balance  horizontal  flow  divergence  amd  conserve  mass.  While  this  vertical  uniformity  assumption 
is  appropriate  in  Cartesiam  coordinates,  it  introduces  spurious  shear  forces  in  a  sphericad  layer.  It 
is  necessary  to  assume  that  the  horizontad  velocity  has  a  radial  dependence.  This  ensures  that 
particles  at  the  top  and  bottom  of  the  layer  that  lie  on  a  Une  through  the  center  of  the  sphere 
have  the  same  amgulair  velocity  about  the  center  of  the  sphere  amd  hence  remaun  in  step  for  any 
horizontal  velocity  field. 


The  appropriate  form  for  an  incompressible  velocity  field  («)  in  a  thin  spherical  layer  is: 


where  A  &  B  are  chosen  to  match  the  conditions: 

w{<f>,6,r  =  a,t)  =  Wb{<f),6,t), 


w{<f),9,r  =  a  +  h,t)  =  Wb{(j),9,t)  +  — . 

at 


This  linear  radial  dependence  of  the  horizontal  components  of  the  velocity  field  is  necessary  to  keep 
the  top  of  a  column  of  fluid  positioned  over  its  bottom  for  all  horizontad  flows.  The  entire  layer 
moves  coherently  as  it  flows  over  the  surface  of  the  sphere  amd  there  is  no  shear  between  fluid  at 
different  radii.  The  verticad  velocity  is  chosen  to  match  the  vertical  motion  of  the  lower  boundary 
amd  the  local  divergence  of  the  horizontal  motion  field.  This  maintains  the  incompressibility  of  the 
overall  flow. 


For  reference: 

w{(l>,9,r,t) 


(r^  -  g^)  (g  +  ^  +  [(r^  +  ha^)  (a  +  -  r^c^]  Wb 

[(a  +  h)^  —  o®] 


(3) 
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In  the  limit  1: 


AT^->rB  n  ,  (r-a\dh 

— -s - y  Wb{4>,9,t)  +  (  , 


h  )  dt' 


similar  to  the  Cartesiain  result  [25]. 
The  leading  error  term  is  : 


-(!)-]■ 


This  vanishes  at  r  =  a  aind  r  =  a  +  h  and  is  bounded  in  that  interval  by: 


The  divergence  of  this  velocity  field  (l)is: 

1  r  d{v'  cos  6)  1  3(a  +  dh  3  [(a  +  —  o^j 

acos9  .  d<t>  ^  d9  .  (o  +  h)^  — dt  ^  (o  +  /»)^  — 

As  this  relationship  does  not  depend  upon  the  vertical  coordinate,  it  is  valid  everywhere  in  the 
fluid  layer. 

In  the  limit  — >  0,  we  recover  (for  an  incompressible  fluid)  a  two-dimensionad  continuity 

equation  : 


dh  1  /  d{hu')  d{hv'  cos0) 

dt  a(X)s9  \  d9 


dh 

|^+V2V«0,  (9) 

(using  the  subscript  2  to  denote  spherical  surface  two  dimensional  vector  calculus  operators)  where 
V  is  the  horizontal  velocity  transport  vector  (see  next  page). 


The  leading  error  term  is 


when  approximating  (7)  by  (9).  In  the  absence  of  vertical  motion  of  the  lower  botmdary,  (9)  is 
aiccurate  to  O  j  . 
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2.2  A  Single  Layer  with  Uniform  Density 

Here  are  the  equations  used  within  a  single  layer.  If  this  is  one  of  a  stack  of  layers,  interactions 
between  layers  axe  mediated  by  the  “sources  -  sinks  ”  terms  in  each  equation.  The  effective  pressure 
terms  ||  will  have  to  be  modified  as  in  NOARL35  [25]  to  reflect  the  possibilties  of  barotropic 
and  baroclinic  modes  of  motion.  Hereafter  we  will  drop  the  prime  to  distinguish  2-D  velocities  from 
3-D  velocities  and  assume  that  (u,  v)  label  horizontal  layer  averaged  velocity  components. 

2.2.1  Continuity 


^  1  cosg) 

dt  ^  a  cos  9  .d4>^  d9 


=  somrces  -  sinks. 


(11) 


2.2.2  Velocity  Transport 

Define  an  angular  deformation  tensor  e  to  have  components: 


11 

(  u  \ 
\cois9) 

0 

<  u  \ 
\c0s9)  ’ 

1 

II 

(12a) 

^(pd  = 

8 

8(f> 

(—) 
\cos9  J 

Q 

+  cos«^ 

(—)■ 
\cos9  J 

II 

(12b) 

where  the  velocities  (u,  v)  are  defined  from  the  transports  (U,  V)  by: 


hu  =  U,  (13a) 

hv  =  V.  (13b) 


Then  the  two  components  of  the  Velocity  TVansport  equation  may  be  written: 
component: 


1 


m 

dt  a  cos  9 


d{Uu)  d{VuCOs9)  .  nrr  •  na 

'  -  -I-  — ^ - -  —  V u  sm0  —  ofi V  sin2d 


d<l>  89 


hg  dh 


d{hcos9e^)  ,  d{hcos^9e<j, 


flcos0  d<j) 

+  sources  -  sinks,  (14a) 


ej  component: 

dV  1 

+ 


8t  a  cos  9 


-|-  hFg  -f- 


Ah 


cos^  9 


d{U  v)  d{v  V  COs9)  .  r^rr  •  nn 

-I-  — — 7^ - -  +  Uu  sin0  -f-  ailU  sm20 

d<f>  89 

8{h  cos  9  e^)  8{h  cos^  9  Ctf,. 

8^  89 


hg  8h 
a  89 


•f  sources  -  sinks. 


(14b) 
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2.2.3  Scalar  transport 


d{hT)  1  r  d{UT)  djVT  cos 6) 

dt  a  cos  6  _  d<f)  ^  dO 


K 

o2  cos^  6 


+  COS0 


)] 


+  sources  -  sinks. 


(15) 


Here  h  is  the  vertical  layer  thickness;  Cl  is  the  rotation  rate  of  the  earth;  g  is  the  acceleration  due 
to  gravity;  and  F$  are  the  components  of  the  apphed  forcing  fields  in  local  {<f>,  0)  components; 
Aff  is  the  eddy  viscosity;  T  is  a  scalar  tracer;  and  k  is  the  coeflicient  of  diff^usivity  for  the  scalar  T. 


2.3  Discussion 


These  forms  are  mathematically  equivalent  to  those  used  in  Shchepetkin  and  O’Brien  (1996) 
[18].  The  form  of  the  dififusion  operator  is  taken  from  Wajsowicz  (1994)  [24],  assuming  that  the 
verticad  diffusivity  and  second  viscosity  coefficients  vanish.  This  form  of  the  friction  term  should 
be  more  accurate  for  layers  with  varying  thickness  than  the  form  used  in  some  earlier  Cartesian 
studies:  AffV^  U ,  in  Hmlbmrt  and  Thompson  (1980)  [8],  although  it  will  require  more  work  to 
calculate. 


Equations  (12a-b)  imply  a  stress-stradn  relationship  of  the  form 


a  =  Anh  cos9  e, 


(16) 


and  a  stress  force  contribution  to  the  transport  equations  of  the  form: 

1 


o  cos  9 


V2 


where  V2'  is  defined  in  Appendix  B  (B5).  of  and  e  are  angular  tensors.  These  require  a  factor 
of  j  to  transform  them  to  physical  tensors  [22].  The  components  of  e  may  be  related  to  the 

rate-of-strain  tensor  S  components  (as  defined  in  [3])  for  the  velocity  field  (1): 


=  o(^^  =  2o5^.  (17) 

Equations  (11)  -  (15)  are  valid  only  in  the  limits  1  and  u,  v  ^  u;.  Typical  terms  ignored 

when  deriving  these  simplified  equations  are  of  the  order  of 
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2.4  Enstrophy  Conserving  Formulation 


An  alternative  formulation  (as  recommended  by  [17],  [1])  uses  the  radial  vorticity  component 
C  defined  as  : 


1  (  dv  d{ucos9) 


a  cos  6 


f  dv 


and  exploits  the  vector  identity 


(u  •  V)  u  =  (V  X  u)  X  i?  +  V 


/u-u 

V~2“ 


This  substitution  yields: 


dU  1  r  dU  dV  d(Vcos0) 

dt  acosy  L  d<p  d(p  d9 


-  (C  +  2f2sin0)  V, 


dV  I  dV  ^  dU  ,  V  tan9Vv]  „ 

H —  +  - - - +  (C  +  2^lsm9)U, 

dt  a  [  d9  d9  cos  9  a  } 

as  an  alternative  left  hand  side  for  equations  (14a)  and  (14b). 


If  u  is  used  instead  of  F  as  the  prognostic  variable,  then  the  left  hand  sides  of  these  equations 
assume  the  simpler  forms  [1]: 


,  dv  h  d  f  +  v^\  ,  .  . 
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2.5  Relocating  the  Poles 

Converging  coordinate  lines  may  cause  difficulties  with  many  numerical  schemes,  particularly 
in  limiting  the  size  of  the  permissible  time  steps.  A  ploy  to  lessen  this  problem  is  to  move  the 

coordinate  singularities  so  that  they  both  lie  over  laind  [19].  If  this  is  done,  the  Coriolis  term  in  the 

momentum  transport  equations  needs  to  reflect  this  loss  of  coincidence  of  the  rotation  axis  and  the 
North  eind  South  Poles  of  the  coordinate  system. 

In  vector  form  the  Coriohs  pseudo-force  in  the  velocity  transport  equations  (14a,  14b)  is  ex¬ 
pressed  as  2  X  U .  Using 

=  ^cos^  sin^n  i  +  sin(^n  sin%  j  -I-  sin^p 

and  re-expressing  in  the  local  spherical  basis  (using  equations  (Ala-c)  from  Appendix  A)  we  find: 

n  =  fi  cos  0n  sin  (0n  —  (p) 

-f  17(sin0n  cos<f>  —  cos^n  sin0  cos{4>q  —  <f>))  eg 
•f  Q  (sin^n  sin^  +  cos^n  cos9  co8{4>u  —  4>) )  ©n 

=  (fl^,  ilg,  Qr) . 

For  U  =  {U,  V,  0)  the  Coriolis  pseudo-force  vector  becomes  2  (— U  flri  U  Clr,  V  —  U  fie) 
or,  in  spherical  surface  coordinates,  (— F  /,  U  f)  where: 

/  =  2fl  (sin^n  sin0  -t-  cos^n  cos9  cos{<f>ii  —  (p)) .  (22) 

Thus,  on  a  shifted  coordinate  spherical  surface,  the  Coriolis  term  becomes  dependent  on  both  the 
local  longitude  and  latitude,  but  retains  its  2fl  {—V  /,  U  f)  form. 
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3.  FINITE  DIFFERENCE  APPROXIMATIONS  TO  THE  EQUATIONS 
3.1  The  grid  on  a  sphere 

Finite  diflFerence  approximations  to  the  equations  set  out  in  Section  2  are  defined  on  a  uniform 
(A<f>,  A0)  grid  covering  most  of  the  earth,  but  commonly  excluding  extreme  latitudes.  The  ’C’ 
grid  of  [13]  defines  the  spatial  staggered  grids  where  the  k,  U  and  V  variables  take  discrete  values. 
Figure  1  sets  out  the  indexing  of  these  variables  on  the  three  interlocked  staggered  meshes. 
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Figure  1 

The  relative  layout  of  the  h,  U  and  V  grids. 

The  dashes  mark  the  boundary  of  the  fluid. 

The  normal  components  of  velocity  vanish  across  this  line. 

A  fourth  grid,  vorticity  (Cu),  can  be  defined  at  the  points  (621+1 ,62i+i),  lying  horizontally  between 

2  2 

the  Vij  points  and  vertically  between  the  Uij  points. 


10 


Moore  and  Wallcraft 


S.1.1  Infinitesmal  vs.  finite  area  elements  and  line  lengths 

An  infinitesmal  difierentiaJ  surface  area  element  on  the  sphere  in  our  coordinates  is 
cos  6  d9d<f).  The  true  area  of  a  finite  mesh  cell  centered  at  {4>,0)  of  width  a  cos6A(f>,  and 
height  aA6  is  2  a?  A<f>  sin  cos  6  .  To  leading  order,  the  proportioned  error  between  the  true 
area  and  the  finite  area  element  a?  cos  9A<f>  A9  is; 

true  area  -  finite  element  A9^  ,  , 

- “  IT- 


More  elaborate  finite  difference  schemes  than  those  set  out  in  the  subsequent  subsections  of 
this  paper  that  intend  to  achieve  higher  accurzicy  will  have  to  incorporate  this  second  order  finite 
area  difference. 

Fortunately,  finite  line  elements  parallel  to  the  (d>,  9)  coordinate  axes  remain  analogs  of  the 
infinitesmal  line  elements;  (a  cos  9  d<j>,  a  d9)  becoming  exactly  (a  cos  9  A<f),  a  A9)  . 

3.1.2  Point  vs.  Area  or  line  averaged  values 

The  variable  hij  is  a  layer  thickness  value  centered  at  the  point  4>i  —  <i>o  +  i  A(f> , 

9j  =  9q  +  j  A9  averaged  over  a  mesh  cell  of  width  acos9jA(l>  and  height  aA9  .  The  variable 
is  an  average  of  a  tramsport  over  a  line  from  the  point  i<l>i,9j)  to  the  point  {(f)i,9j+i)  . 
Similarly,  is  defined  as  the  average  of  the  transport  across  the  line  from  {<f>i,9j)  to  the 

point  {(f>i+i,9j) .  In  our  spherical  coordinate  system  these  averages  differ  from  spot  values  centered 
on  {(f>i,9j)  by  terms  of  second  order  and  higher  in  A9,A(p  times  and  9  derivatives  of  the 
field. 


If  we  consider  in  detail  the  relationship  between  a  spot  vaJue  of  a  field  h{d>o,9o)  and  the 
average  value  of  this  field  over  the  finite  area  element  of  size  2o^  Aip  sin  cos  9  centered  at 

the  point  {4>q,  9o)  {called  h  )  we  find: 


d^h 

(A0)2 

'd?h 

«  ..  dh 

‘ 

+  — — 
-  24 

d9^ 

-2  tan^o  ^ 

^0 

00  , 

A4?A9^,A9^) 


+ 


(24) 


A  similar  second  order  dependence  on  the  mesh  spacing  exists  between  the  average  value  of 
U  between  the  points  {<pi,9j)  and  {<f>i^9j+i)  and  its  spot  value  U{4>i,9j^i).  Likewise  for  V. 

Again  these  differences  between  spot  and  average  values  must  be  reflected  in  more  accurate 
difference  formulae. 
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3.2  Second  Order  Finite  Difference  Equations 

For  convenience  in  notation,  let  ns  follow  other  authors  in  ocean  modelling  [11]  and  define 
centered  averaging  and  difference  operators: 

Smz  {W{Z))  =  -w(z-  /(mA),  (25a) 

=^w(z-\-^"j  +W(z-  2.  (25b) 


The  centered  time  t  applies  to  all  fields  unless  otherwise  noted.  Second  order  accmrate  finite 
difference  approximations  to  the  equations  set  out  in  Section  2.2  take  the  following  forms: 

3.2.1  Conservation  of  Mass 

Equation  (11): 


S2th  = 


-1 

acosO 


S4,U  +  6g{V  cos 9) 


+  sources  —  sinks. 


(26) 


3.2.2  Velocities 


The  velocities  (u,  v)  are  defined  from  the  transports  (U,  V)  by: 

Vi. 


■■ 


Vt,j  -  .  . 


(27) 


3.2.3  Deformation  Tensor  Components 
Equations  (12a)  and  (12b): 


,  (28a) 

=  e,,  =  J*(j^)  +  cosWs(^),  (28b) 

This  defines  the  diagonal  elements  of  the  deformation  tensor  e^  and  egg  on  the  hi^  grid  and 

the  off-diagonal  elements  e^^g  2md  eg^  on  the  Ctj  grid-  This  positions  these  values  exactly  where 

they  are  needed  by  the  difference  operators  to  calculate  the  viscous  stress  contribution  to  the  U 
and  V  advection  equations. 
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3.2.4  Velocity  transport 

Equations  (14a)  and  (14b)  are  approximated  by: 

^2tU  =  — +  6g(v‘^rf  cos  9^  —  sin^  (u^/i)^u  —  a  Cl  sin  20  (u^ /i)*^ 

T<t>  _  1/  r  — 1 1  -  At 

--^^64,h+  ^ ^  64h  cos  9  e^)  +  SeChT  cos"^  9  e^)  ,  (29a) 

a  COS  6  ^  a'cos^0  L  \ 


62tV  =  S^(Tf¥)  +  6e(y^lf  cos9)+  {sin9Uu^h)  +  of!  (sin  29trt>hf 

-  ^Seh+  jfFg  +  ^2cos^0  k(^cos0efl^)  -  6e[hcos^ 9e^)\  .  (29b) 

3.2.5  Scalar  tracers 

Scalar  fields  such  as  temperature,  salinity  and  tracer  concentrations  (T,  5,  Ci)  are  defined  on 
the  hij  grid.  Equation  (15)  may  be  approximated  in  a  manner  consistent  with  the  previous  finite 
difference  equations  by: 

ht{hT)  =  WiUT^)  +  SeiVT^  COS  5)1 

a  cos  u  I 

+  S^Tj  +  cos 9 6e(hcos 9° 6gT)  (30) 

3.2.6  Discussion 

The  horizontal  diffusion  terms  in  the  three  finite  difference  equations  given  above  are  lagged  in 
time  one  time  step  for  computational  stability  [16].  This  introduces  first  order  timestepping  errors. 
In  practice,  small  values  of  u  and  k  are  used,  in  fact  the  smallest  that  allow  the  code  to  be  stable. 
Hence  a  small  constant  proceeds  this  first  order  timestep  error  amd  its  overall  effect  is  hoped  to  be 
small. 

In  the  momentum  equations  for  time  stepping  U  Sc  V,  the  advection  terms  have  a  conserva¬ 
tive  fiux  contribution  and  separate  curvature  terms  (—V  u  tain0,  U  u  tan0).  These  curvature  terms 
redistribute  the  momentum,  but  if  the  equation  {u  x  (14a)  -f  u  x  (14b))  is  formed  to  pro¬ 
duce  a  prognostic  equation  for  a  kinetic  energy  {uU  -I-  v  F),  it  is  seen  that  these  curvature  terms 
make  no  contribution  to  the  kinetic  energy  evolution,  just  as  the  Coriolis  terms  also  cancel  in  this 
combination. 
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Holland  and  Lin  [5]  devised  finite  difierence  forms  for  the  Coriolis  terms  in  each 
equation  {—{v^fh)  ,{u’f’  fh)  )  that  preserve  this  conservation  of  kinetic  energy  when 
f  h)  is  summed  over  a  domain  where  the  normal  velocity  component  van¬ 


ishes  at  the  boundary.  This  may  be  confirmed  numerically  for  arbitrary  (/,  h,  U,  V)  fields  satisfying 
the  boundary  conditions.  When  the  axis  of  rotation  goes  through  the  coordinate  singularities 
(North  &  South  Poles),  this  finite  difierence  approximation  reduces  to  the  form  set  out  above. 
The  more  general  forms  given  in  this  paragraph  should  be  used  when  the  poles  are  moved.  Other 
suggested  finite  difierence  approximations  for  the  Coriolis  terms  [18]  do  not  conserve  kinetic  energy. 


This  energy  conservative  finite  difierence  approximation  to  the  Coriolis  terms  is  chosen  as  a 
guide  to  approximating  the  curvature  term  in  each  equation. 


The  scalar  tracer  finite  difierence  equation  preserves  locally  the  property  of  no  change  in  the 
scalar  field  when  it  is  constant,  no  matter  what  the  velocity  field.  Any  local  divergence  in  the  17,  V 
field  is  reflected  only  in  a  change  of  layer  thickness. 


3.3  Diagnostic  Quantities  and  Energetics 

3.3.1  Vorticity 


The  normal  component  of  vorticity  is  defined  on  a  spherical  surface  grid  [4]  as: 

.  _  1  f  ^  d{u  cos  0)  "1 

a  cos  6  \  d(l>  do  f  ' 


The  finite  difierence  analogue  of  this  equation  is: 

Cm  =  -  <^e(«cose)}. 


(32) 


It  may  be  used  as  a  diagnostic  variable. 

3.3.2  Kinetic  Energy 

If  the  kinetic  energy  field  is  defined  as  Uij  Uij  +  Vjj  Vij,  the  two  products  are  defined  on 
separate  grids.  One  or  both  need  to  be  extrapolated  onto  a  common  grid  to  define  a  kinetic  energy 
field,  although  the  separate  grids  may  be  used  for  the  purposes  of  integrating  kinetic  energy  over 
the  entire  domain.  If  the  kinetic  energy  is  to  be  combined  or  compared  with  the  potential  energy, 
then  the  definition  on  the  /ijj  grid  would  be  appropriate,  viz.\ 

{KE)ij  =  Uij  UiJ^  +  Vij  Vij  . 


(33) 
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4.  TRUNCATION  ERROR  AND  ACCURACY 

As  the  operators  (25a, 25b)  used  in  constructing  the  finite  difference  approximations  to  the 
equations  of  fluid  motion  are  centered,  we  expect  them  to  represent  these  differential  equations 
with  error  terms  proportional  to  even  powers  of  the  mesh  spsicings  (A<^,  A0)  and  their  products. 
A  full  truncation  error  analysis  of  the  finite  difference  equations  in  section  3.2  will  reveal  this 
expected  behavior. 

It  is  useful  to  attempt  to  quantify  these  error  expressions  by  assuming  a  surface  flow  of  a 
simple  type  and  seeing  how  these  finite  difference  equations  reproduce  the  properties  of  the  analjrtic 
equations  for  this  flow.  The  explicit  form  taken  by  the  truncation  errors  in  this  case  may  be  a  useful 
guide  as  to  the  choice  of  mesh  spacing  and  to  the  degree  the  distortion  of  the  mesh  near  the  poles 
degrades  the  accuracy  of  these  equations.  One  such  simple  flow  is  that  of  a  rigid  shell  in  unform 
rotation  about  some  pole  centered  at  the  point  0^) . 

4.1  Analytic 

The  surface  velocity  field  (u,  u]  of  a  thin  spherical  shell  rotating  firom  west  to  east  about  an 
axis  passing  through  the  point  {<1>r.  Or)  with  angular  velocity  lj  is: 

Uii  =  a  a;  [cos  (<^  —  (f>R)  cos  sin  ^  —  sin^ijcos^,  sin{<pR  —  4>)  cos0/j].  (34) 


If  ua  is  used  to  prescribe  U  and  V  in  the  continuity  equation  (11),  it  can  be  shown  that  the 
expression  in  the  square  brackets  vanishes  identically.  If  Uij  is  inserted  into  the  expressions  for  the 
viscous  stress  tensor  components  then  it  be  shown  that  they  aJl  vanish  everywhere,  as  expected  for 
such  a  flow. 

If  these  forms  are  used  for  U,  V,  u  and  v  in  the  bilinear  advection  terms  for  the  velocity 
transport  time  evolution,  the  terms  that  arise  can  be  balanced  exactly  by  the  pseudo-Coriolis  force 
for  an  axis  of  rotation  moved  to  ( <f>R,  Sr)  and  a  rotation  rate  of  —  w .  This  is  as  expected  when 
the  apparent  forces  arising  &om  a  rotation  of  the  inertial  axes  is  cancelled  by  an  apparent  flow  field 
in  the  moving  frame  exactly  opposite  to  the  rotation  field. 

These  analytic  results  will  be  useful  for  investigating  the  accuracy  of  finite  difference  represen¬ 
tations  of  the  governing  equations  and  alternative  formulations. 
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4.2  Numerical 


If  the  above  velocity  field  is  projected  onto  a  grid  of  size  (A^,  Aff)  and  these  discrete  values 
are  used  in  the  finite  difference  approximations  for  the  governing  equations  we  can  look  at  the 
differences  between  the  discrete  results  and  the  analytic  results  in  our  several  cases. 


4-^-J  Continuity 


The  right  hand  side  of  equation  (26)  evaluates  to: 


cos  Or  sin  <f)  tan  9 


(35) 


Thus  if  =  A9 ,  the  terms  inside  the  square  bracket  cancel  and  our  finite  difference  scheme  is 
exact  for  this  model  flow.  The  difference  between  A</>  and  A0  determines  the  actual  truncation 
error  for  this  scheme  and  grid.  To  the  lowest  order,  this  error  takes  the  form: 


62t  h  «  cos  9r  sin  4)  tan  6 


(M! 

24 


^  +  O((A#-(A0)^)  +  ... 


4.2.S  Advection  terms 

K  the  general  rigid  rotation  velocity  field  (34)  is  substituted  into  the  advection  part  of  the 
trainsport  equations,  very  complicated  expressions  result.  The  general  level  of  accuracy  of  the  finite 
difference  approximations  here  can  be  illustrated  by  considering  two  limits  of  the  rotation  velocity 
field;  (i)  axis  of  rotation  at  the  poles  and  (ii)  axis  of  rotation  at  the  equator. 

Let  the  velocity  field  for  case  (i)  be  «  =  (—  cos  0,  0 ) .  This  velocity  field  in  the  co-rotating 
firame  exactly  recovers  the  inertial  frame  and  there  are  no  real  forces  in  the  zidvection  terms.  By 
inspection  we  can  see  that  the  finite  difference  approximations  to  the  advection  terms  in  the  U 
transport  equation  (29a)  evaluate  to  zero  for  any  choice  of  A4>  and  A6 .  In  the  V  tranport 
approximation  (29b)  the  curvature  term  is  cancelled  exactly  for  any  mesh  by  the  Coriolis  force 
finite  difference  approximation. 


When  we  move  the  axis  of  rotation  to  the  equator  at  =  0 ,  the  velocity  field  takes  the 
form:  u  =  ( —  cos(^  sin0,  sin<^) .  The  analytic  contribution  to  the  U  transport  advection  is 
sin<^  cos  <!>  cos^  6 .  The  finite  difference  approximation  evciluates  to: 


sin  4>  cos  <f)  cos^  6 


,  3  +  \Ais.v?9  .  ,2 

1 - 24 - 


7  cos  20 
24  cos^  9 


A9^ 


+ 


(36) 


Both  of  the  coefficients  of  the  second  order  terms  axe  well  behaved  between  70®  N  and  70®  S ,  being 
bounded  by  3.3  for  the  first  coefficient  and  -1.25  for  the  second  in  this  range.  The  truncation  error 
properties  of  the  finite  difference  approximation  to  the  V  transport  equation  (29b)  are  similar. 


16 


Moore  and  Wallcraft 


4.2.3  Deformation  tensor 

iProm  the  components  (25a-b)  for  the  deformation  tensor,  we  recover  truncation  error  expres¬ 
sions  of  the  form 

e^(uii)  «  cos 9r  sin<j)  tamO  — I-  ^5  -f  6  tan^(^)^  — h 


e,^(uit)  «  cos 6r  cos (i>  secO 

Although  the  factors  involving  6  tern  ^6  can  get  leirge  at  high  latitudes,  these  approximations  are 
seen  to  retaun  very  good  second  order  accuracy  over  much  of  the  globe. 

When  these  expressions  are  used  in  calculating  the  stress  contribution  to  the  velocity  trans¬ 
port  equations  (26a-b),  the  resultant  truncation  errors  retain  their  ( (A^)^,  {A6)^ )  character  and 
magnitude. 
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Appendix  A 


THE  CONNECTION  BETWEEN  CARTESIAN  AND  SPHERICAL  BASIS 

VECTORS 


It  is  useful  to  be  able  to  re-express  vectors  defined  in  an  absolute  (i,  j,  k)  Cartesian  basis  in 
terms  of  the  local  spherical  basis  ( e,^,  eg,  Cr)  at  the  point  {<f>,  6,  r)  eind  vice  versa.  The  nine 
components  of  the  matrix  that  transforms  a  vector  representation  from  one  basis  to  the  other  can 
be  expressed  in  either  the  local  spherical  coordinates  {<f>,  6,  r)  or  in  the  local  cartesian  coordinates 
(x,  y,  z).  We  will  assume  for  convenience  that  the  origins  of  the  two  coordinate  systems  coincide 
and  that  the  two  spherical  coordinate  singularities  (the  North  and  South  Poles)  lie  on  the  0-axis. 

A.l  Spherical  coordinate  representations 

To  move  from  a  Cartesian  basis  to  a  spherical  basis  we  use  the  relations: 


=  —  sm<f) 

-t 

i  + 

C0s4> 

J 

(Ala) 

eg 

=  —cos4>sm6 

• 

1  — 

sincf)  sind 

j  -b  COS  9  ic, 

(Alb) 

St 

=  cos  <f>  cos  0 

i  + 

sin<^  cos  9 

j  +  sin^  ic. 

(Ale) 

To  move  from  a  representation  in  the  spherical  basis  to  a  representation  in  the  Cartesian  basis 
we  use  the  relations: 

i  =  —  sin<^  e^  —  cos<^sin0  e^  -1-  cosi^cos^  Sr 

j  =  cos</>  —  sin^sin0  eg  -I-  sin^cos0  Sr 

k  =  COS0  S^  -I-  sin^  Sr. 

If  these  transformation  rules  are  expressed  as  Vj  =  T  vc,  v,  is  the  3  components  of  v  expressed 
in  the  local  spherical  basis  set  {  e^,  eg,  Sr),  vc  is  that  same  vector  expressed  in  the  Cartesian 

basis  set  (i,  j,  k^  and  T  is  a  three  by  three  matrix.  We  note  that  ^T^  =  ^T^  .  The  matrix 

T  is  orthogonal  as  its  transpose  is  its  inverse! 


(A2a) 

(A2b) 

(A2c) 
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A.2  The  Spherical-Cartesian  transformations  in  Cartesian  Coordinates 


P  = 


To  express  T  purely  in  terms  of  local  Cartesian  variables,  we  define  r  =  {x^  +  + 

(x^  +  Then 


1 

2 

y 


X 

p 

rp 


\ 


X 

r 


k 

r 


) 


(A3) 


In  this  form  it  is  easy  to  demonstrate  that  the  basis  (e^,  eg,  Cr)  defined  by  (A1  a-c)  has  the 
(required  !)  properties: 


e^  = 

1, 

e^  •  eg 

=  0, 

e0*  er 

=  0, 

X 

4  = 

®r) 

ee  •  eg 

=  1, 

•  Cp 

=  0, 

(A4) 

Br  X 

e^  = 

ee 

Cr  X  eg 

60, 

Bp  •  Bp 

=  1. 

Appendix  B 

SPHERICAL  SHELL  VECTOR  CALCULUS 


B.l  Definition  of  Surface  Differential  Operators  on  a  Sphere 
B.l.a  Basic  Operators 

Let  $  be  any  scalar  quantity,  and  let  it  be  a  function  of  (0, 6)  only, 

(mz.  ^  =  0). 

Let  P  be  a  vector  quantity  such  as  velocity  v,  or  mass  transport  V.  F  is  assumed  to  lie  in  the 
local  tangent  plane  and  to  describe  a  field  ever3rwhere  nearly  parallel  to  the  spherical  layer  it  is 
embedded  in.  Quantities  advected  by  v  {or  V)  are  assumed  to  maintain  vertical  coherence  across 
a  fluid  layer  (c.  f.  [14]  p  63).  That  is,  two  particles  lying  on  the  same  radial  vector  at  the  top  and 
bottom  of  a  fluid  layer  respectively  remain  on  a  single  radial  vector  as  they  are  advected  by  the 
velocity  field  of  that  layer.  This  implies  that  u  is  obliged  to  assume  the  form: 

u(<^,  e,r)=  (-  u'{4>,  9),  -  v'{(f),  9),  w{<f),  0,  r)  V  (Bla) 

\  a  a  J 

where  w{4>,  9,  r)  takes  the  form  necessary  to  maintain  zero  divergence  for  the  velocity  and  to 

match  the  motions  of  the  upper  and  lower  boundaries  of  the  layer,  in  what  follows,  let  us  consider 
a  general  vector  F  is  expressed  in  coordinate  form  as: 

F  =  '^F^{<l>,9),Fr{<f>,9,r)'^  .  (Bib) 

In  many  ocean  model  contexts,  the  radial  component  of  the  velocity  vector  is  ignored.  For  these 
purposes  it  is  useful  to  consider  just  the  horizontal*  components  of  velocity.  Let  us  define  U2  to 
be: 


U2{(t>,0,r)  =  {u’{<f>,9),  v'{4>,9),  0),  (B2) 

and  similarly  F^  to  be  any  locally  horizontal  vector  field. 

The  standard  invariant  operators  of  vector  calculus  are  now  set  out  in  explicit  coordinate  form 
using  our  spherical  surface  variables  {cp,  9 )  and  assuming  a  functional  form  for  vectors  defined  as 
in  (Bib)  or  (B2)  and  for  scalars  dependent  upon  {(f>,9)  only.^ 


^c.f.  Batchelor  (1970),  pp  598-601  or  Gill  (1982),  pp  92-93 
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B.l.b  Gradient  (V2^) 


*  r'Se‘>- 


(B3) 


This  operator  retains  an  r  dependence  throughout  a  spherical  shell  and  will  have  to  be  vertically 
averaged  with  care. 

B.l.c  Divergence  (V  •  F) 


V-F  =  1  djFjcosd)  1  djr^Fr) 

a  COS  6  d<f>  a  cos  6  dO  dr 


(B4) 


The  first  two  terms  are  independent  of  r.  This  expression  may  be  used  to  define  an  r  component 
to  a  vector  field  to  madce  it  divergence  firee.  We  will  define  V2*  to  be: 


V2  •  F2  = 


1  ^  1  5(F^cos0) 

a  cos  6  d<t>  acosO  89 


(B5) 


B.l.d 


Curl  (V  X  F) 


V  X  F 


=  ri  V  _ l_  ^ 

Vr  5^  a  )  y  a  rcosO  d4> 


+ 


1  f  ^  _  g(F^cosg)  1  ^ 
acosd  \  d<i)  89  j 


(B6) 


-J 

When  F2  is  the  velocity  field,  the  Cr  is  the  normal  component  of  the  spherical  surface  flow  vorticity, 
commonly  labeled  Note  that  there  are  both  and  e$  components  of  the  full  vorticity  vector 
for  all  non-zero  horizontal  velocity  fields,  even  in  the  absence  of  any  vertical  velocity. 


B.l.e  Scalar  Laplacian 


(vi$) 


vH  = 


^  1  ^ 


r^cos^9  84^  cos  9  89  \  89 


cos  . 


(B7) 


^Pedlosky  (1987),  p.  55,  eqn.  (2.9.20)  is  incorrect.  Consider  Stokeses  theorem  : 

{V  X  v) '  ds  =  i?  •  for  the  infinitesmal  surface  area  element  A<p  cos  ^  . 
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B.l.f  Vector  Laplacian 

1 


(V2/2) 


= 


ra  cos^  6 
1 

ra  cos^  0 
2 

ra  cos  9 


!(• 


ei^FL 


^  a(F;cosg)1 

d<f>  do 


dF^ 


cos9-^  +  cos29FL  “  2sin0 

0(7  *  ^ 


dFL 


9  /  dF^  \ 

Y  +  cos9—  (cos9-^J  +  cos29Fff  +  2sin0 


d<p 

dF^- 


d(t> 


(B8) 


We  recognize  the  er  component  as  —  ^  V2  •  F2. 

B.2  Higher  Vector  Derivatives  of  the  velocity  field 


B.2.a  Momentum  Advection 

(U2  •  V)  «2  = 


It 


r 

+  - 


dv!  v'  du'  tan0  ,  , 

a  Lacos0  d(f)  a  89  a 

dv'  ^  v'  dv'  ^  tan6  ,  / 
a  Locos0  d<f)  a  89  a 

^u'u'  +  v'v' 

■  er 


u 


a 


(B9) 


B.2.b  curl  curl 

V  X  (V  X  U2)  =  - 


+ 


ra  cos2  9 
1 

ra  cos2  9 
2 

r a  cos  9 


9  /  9‘ix^  \  9t^^ 

COS0—  (cos^-;^  1  +  cos20u'  —  sin0 

V  89  J 


89 


8(f) 


—  cos  9 


^v' 


^  +  2cos^«t.'-  ^(cose^) 


cf^v 
8(f>^ 

8u'  8(v'cos9) 


ee 


er. 


8^89 


(BIO) 


Note:  V  X  ( V  X  «2)  /  —  V2  U2  because  V2  •  tt2  #  0 

B.2.C  grad  div 

1 


V  (V  •t?2)  =  + 
+ 


ra  cos2  9 
1 


8\'  8  (  ^ 
8(p  ^  89 


ra 

+  0-  er 


{—[ 

1 

.89  ' 

Vcos0  \ 

8v 

T4> 

8{v'  cos  9) 


D] 


^  '  d9 


(Bll) 
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B.2.d  div  grad 


(V-V)U2  = 


1  r  a  ^  (  a  - 

+  0  *  er 


(B12) 


B.3  The  rate-of-strain  tensor 


Assuming  w  is  not  zero: 


S(f^ 

£e<t>  ^69  ^Or 

^t9  ^rr 


_  du; 


(B13) 


(B14) 


2r  Oe' 


(B15) 


Assuming  w  is  zero: 


c  _c  1 

-  2rcoS^  d^’ 


(B16) 


^99  = 


a  d0’ 


(B17) 


= 


1  dv!  v'  tan  9 
acosO  d<l>  a  ’ 


(B18) 


1  du' 
2ocos0  d^’ 


(B19) 


